#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#include "../src/em.c"

typedef struct _EmTestDatum EmTestDatum;
struct _EmTestDatum 
{ 
  struct
    {
      float    loga[6 * 4];
      float    logbeta;
      float    c[6 * 5];
      float    d[6 * 5];
      float    pz[5];
      float    logpriorz[5];
      int      clamp;
      Rating   ratings[8];
      float    logbetamean;
      float    logbetastddev;
    }                   inputs;

  struct
    {
      float    newq;
      float    logbeta;
      float    pz[5];
    }                   desired_output;
};

static void
test_em (void)
{
  EmTestDatum data[] = 
    {{{{0.753458, 0.0700157, 0.849855, 0.265468, 0.127633, 0.188873,    0.434107, 0.59523, 0.196637, 0.783694, 0.455691, 0.392038,    0.736491, 0.180422, 0.887738, 0.754236, 0.692066, 0.66938,    0.074632, 0.453754, 0.235292, 0.985654, 0.225451, 0.343979},   0.481835, {0.915638, 0.375595, 0.078511, 0.354202, 0.726765,    0.941488, 0.483281, 0.157565, 0.943072, 0.485797, 0.0912432,    0.421074, 0.76265, 0.598058, 0.337007, 0.729008, 0.0932695,    0.523427, 0.883253, 0.493716, 0.107615, 0.297976, 0.539275,    0.011881, 0.191977, 0.922381, 0.460764, 0.657679, 0.465212,    0.980893}, {0.977482, 0.500114, 0.52214, 0.495096, 0.886239,    0.0790406, 0.75949, 0.897037, 0.549232, 0.350033, 0.666221,    0.373611, 0.665979, 0.856317, 0.558605, 0.075635, 0.126704,    0.844436, 0.366628, 0.153254, 0.66594, 0.186757, 0.901416,    0.172362, 0.688458, 0.686643, 0.379276, 0.677266, 0.802219,    0.607602}, {0.619786, 0.780228, 0.252987, 0.257569,    0.953566}, {0.406617, 0.587008, 0.401252, 0.39496,    0.330982}, -1, {{3, 2}, {2, 3}, {2, 3}, {1, 3}, {0, 3}, {0,     0}, {3, 0}, {3, 3}}, 0.556816,   0.0283322}, {-17.278206721077845088,   0.54124214449138466713, {1.7463e-11, 5.72578e-13,    5.93502e-12, 1., 4.65651e-16}}}, {{{0.177728, 0.794364,    0.370059, 0.126916, 0.0053663, 0.105906, 0.683416, 0.747639,    0.328101, 0.303687, 0.0758142, 0.127853, 0.547873, 0.0507004,    0.818245, 0.174288, 0.141255, 0.463692, 0.416993, 0.779327,    0.810273, 0.00338782, 0.860177, 0.750995},   0.632545, {0.209024, 0.490118, 0.624079, 0.627179, 0.103117,    0.806701, 0.87644, 0.299078, 0.79943, 0.730887, 0.748587,    0.751206, 0.74873, 0.912642, 0.574299, 0.60995, 0.285037,    0.495649, 0.794972, 0.799677, 0.28165, 0.635473, 0.0439768,    0.167132, 0.072626, 0.145355, 0.419897, 0.539952, 0.969509,    0.338654}, {0.543457, 0.240874, 0.170078, 0.607767, 0.79487,    0.489668, 0.421349, 0.695125, 0.220571, 0.879718, 0.136311,    0.199476, 0.425599, 0.0800407, 0.854662, 0.564003, 0.381622,    0.912909, 0.782036, 0.418647, 0.961725, 0.372957, 0.812527,    0.0799933, 0.418268, 0.132083, 0.642449, 0.472226, 0.623397,    0.642415}, {0.2211, 0.777101, 0.402826, 0.762697,    0.0847886}, {0.577626, 0.977227, 0.682656, 0.230127, 0.0136232},   0, {{3, 2}, {0, 1}, {3, 3}, {3, 3}, {2, 3}, {2, 0}, {1, 1}, {3,     1}}, 0.769747,   0.448091}, {-21.657043972877618871, -0.65889583406204366947,{0.999875, 0.0000562685, 4.1909e-6, 0.0000644459,    4.89143e-7}}}, {{{0.594976, 0.63388, 0.39679, 0.635564,    0.514983, 0.215612, 0.264707, 0.993115, 0.0427562, 0.592214,    0.622293, 0.772015, 0.265655, 0.189388, 0.859596, 0.687226,    0.688029, 0.212161, 0.176939, 0.457099, 0.674406, 0.616556,    0.407192, 0.0090084},   0.0794299, {0.982677, 0.0104016, 0.373445, 0.564447, 0.767065,    0.745694, 0.38033, 0.521691, 0.174851, 0.123401, 0.608315,    0.256036, 0.985463, 0.263806, 0.921088, 0.568007, 0.773302,    0.0868666, 0.463989, 0.893602, 0.156745, 0.679675, 0.454981,    0.814172, 0.174068, 0.669273, 0.0815359, 0.249724, 0.407003,    0.923579}, {0.701206, 0.728033, 0.232153, 0.800177, 0.0928915,    0.471997, 0.24669, 0.536371, 0.171803, 0.90399, 0.473389,    0.449505, 0.707814, 0.0103878, 0.316643, 0.76983, 0.252834,    0.196216, 0.142575, 0.100557, 0.171298, 0.946492, 0.735572,    0.176978, 0.470092, 0.218458, 0.503419, 0.376801, 0.3772,    0.746461}, {0.256729, 0.84043, 0.205397, 0.842472,    0.78334}, {0.390925, 0.497583, 0.832084, 0.466697,    0.621095}, -1, {{1, 2}, {3, 1}, {0, 1}, {3, 3}, {3, 0}, {1,     3}, {3, 1}, {0, 4}}, 0.635868,   0.324122}, {-18.152851439069185165, -0.21394010038925245653,{0.0003936, 0.969309, 0.005021, 0.0246537,    0.000622918}}}, {{{0.520537, 0.0734509, 0.689376, 0.58855,    0.343559, 0.603359, 0.470918, 0.085131, 0.966758, 0.226159,    0.724457, 0.828402, 0.126328, 0.0207623, 0.881985, 0.0450615,    0.735403, 0.52318, 0.0499015, 0.578364, 0.114309, 0.278431,    0.414034, 0.254243},   0.593771, {0.20498, 0.724657, 0.665692, 0.250212, 0.601621,    0.253739, 0.580561, 0.283455, 0.375461, 0.529282, 0.752159,    0.157127, 0.354699, 0.647297, 0.707098, 0.421723, 0.831519,    0.597396, 0.128733, 0.307415, 0.553089, 0.183362, 0.874491,    0.713643, 0.348109, 0.458704, 0.208799, 0.463431, 0.746488,    0.204965}, {0.628237, 0.179976, 0.371027, 0.675683, 0.876078,    0.0228494, 0.0163277, 0.0283856, 0.16898, 0.601126, 0.184808,    0.43099, 0.0402465, 0.293711, 0.63172, 0.247628, 0.165756,    0.580068, 0.283611, 0.788924, 0.956957, 0.116637, 0.537123,    0.583959, 0.32872, 0.936661, 0.166096, 0.908276, 0.452642,    0.913812}, {0.149768, 0.87989, 0.283662, 0.312686,    0.96496}, {0.4489, 0.243416, 0.0189746, 0.33324,    0.201272}, -1, {{1, 4}, {3, 0}, {1, 2}, {1, 3}, {1, 4}, {3,     3}, {1, 1}, {3, 3}}, 0.438906,   0.0496294}, {-26.338502984496192539,   0.37391154603290714264, {0.0000773747, 0.000613301, 0.0000310435,    0.251524, 0.747754}}}, {{{0.412348, 0.120703, 0.322269, 0.512507,    0.82839, 0.791983, 0.385607, 0.346411, 0.920114, 0.339341,    0.471795, 0.196643, 0.0402231, 0.055679, 0.159109, 0.231683,    0.591323, 0.812263, 0.140135, 0.898442, 0.390051, 0.734604,    0.701228, 0.848813},   0.977702, {0.613901, 0.378959, 0.336306, 0.149313, 0.821918,    0.993352, 0.989896, 0.229199, 0.482577, 0.521557, 0.793253,    0.188976, 0.426898, 0.362448, 0.56157, 0.597653, 0.614634,    0.222313, 0.663128, 0.207603, 0.880031, 0.521085, 0.814315,    0.2299, 0.26613, 0.142125, 0.478009, 0.0805877, 0.444212,    0.148773}, {0.488113, 0.851389, 0.961636, 0.627216, 0.69486,    0.662412, 0.534738, 0.264769, 0.133289, 0.064759, 0.920104,    0.042456, 0.470161, 0.857156, 0.0400732, 0.521371, 0.655846,    0.627256, 0.773943, 0.379246, 0.177838, 0.546668, 0.329731,    0.230473, 0.689725, 0.69528, 0.368095, 0.603256, 0.994865,    0.0328671}, {0.833357, 0.338487, 0.861576, 0.968108,    0.913254}, {0.296031, 0.391415, 0.110952, 0.87318, 0.77466},   4, {{0, 4}, {3, 1}, {3, 1}, {0, 2}, {1, 0}, {3, 0}, {2, 2}, {3,     4}}, 0.483696, 0.0992371}, {-30.541057889589055095,   0.22045907326978127443, {1.00019e-7, 1.38313e-6,    6.18322e-6, 5.87372e-10, 0.999992}}}, {{{0.395414, 0.55773,    0.937028, 0.769506, 0.164941, 0.868005, 0.241749, 0.401411,    0.561685, 0.87314, 0.208881, 0.568053, 0.223197, 0.0115643,    0.240773, 0.6548, 0.927166, 0.62015, 0.129822, 0.78162, 0.152506,    0.884581, 0.646125, 0.682382},   0.757092, {0.326851, 0.709097, 0.912876, 0.592151, 0.458845,    0.467349, 0.511466, 0.0304663, 0.585705, 0.258467, 0.943412,    0.807269, 0.574141, 0.0176941, 0.288612, 0.880103, 0.953991,    0.887873, 0.506993, 0.727597, 0.0694099, 0.241747, 0.82461,    0.970505, 0.742559, 0.53265, 0.911734, 0.378354, 0.283714,    0.0653008}, {0.400268, 0.347888, 0.698009, 0.806833, 0.456856,    0.540619, 0.123868, 0.789139, 0.168244, 0.660516, 0.169877,    0.901267, 0.661251, 0.932919, 0.100467, 0.659519, 0.83664,    0.962414, 0.357908, 0.12687, 0.924906, 0.584059, 0.074194,    0.0615689, 0.524638, 0.236172, 0.376185, 0.254736, 0.0677821,    0.695553}, {0.252317, 0.465596, 0.899538, 0.0350371,    0.0824402}, {0.56433, 0.238288, 0.102118, 0.981973,    0.904811}, -1, {{3, 4}, {3, 3}, {2, 3}, {3, 1}, {2, 1}, {2,     2}, {3, 0}, {2, 4}}, 0.139705,   0.624065}, {-14.807898600100085828, -1.0168432053591173221,{0.0688046, 0.150795, 0.0293288, 0.351907, 0.399165}}}, {{{0.777941,    0.476741, 0.555645, 0.549872, 0.716372, 0.952103, 0.319474,    0.173686, 0.461636, 0.884321, 0.623921, 0.921369, 0.99604,    0.984783, 0.588884, 0.838929, 0.43171, 0.746495, 0.486765,    0.856956, 0.5269, 0.344847, 0.34706, 0.23289},   0.748959, {0.868106, 0.791415, 0.683019, 0.0325867, 0.916003,    0.471941, 0.509332, 0.57095, 0.0316818, 0.84802, 0.587963,    0.57491, 0.0468992, 0.259137, 0.749034, 0.1432, 0.300404,    0.772371, 0.892078, 0.6163, 0.955557, 0.425311, 0.659188,    0.867342, 0.0874511, 0.633896, 0.976169, 0.834755, 0.171448,    0.161955}, {0.466837, 0.263805, 0.139766, 0.313935, 0.878874,    0.688895, 0.0928672, 0.054798, 0.129839, 0.545695, 0.792463,    0.282427, 0.237761, 0.929394, 0.836906, 0.857116, 0.578573,    0.0620527, 0.749455, 0.22322, 0.602404, 0.227298, 0.578006,    0.0612647, 0.135567, 0.963493, 0.43824, 0.74733, 0.256694,    0.274598}, {0.345373, 0.692532, 0.126854, 0.728903,    0.55291}, {0.410105, 0.889093, 0.799509, 0.716004,    0.55299}, -1, {{3, 4}, {1, 0}, {0, 3}, {3, 0}, {1, 1}, {0, 3}, {0,      0}, {0, 0}}, 0.737456,   0.96655}, {-13.541660618889708825, -0.82873067799634964525,{0.974584, 0.00138786, 0.000231903, 0.0225285,    0.00126811}}}, {{{0.32977, 0.708116, 0.510159, 0.388543, 0.268505,     0.572549, 0.546666, 0.950303, 0.521175, 0.315855, 0.272068,    0.604931, 0.828643, 0.189, 0.543164, 0.0520205, 0.418538,    0.299907, 0.743655, 0.336016, 0.865548, 0.989387, 0.00619887,    0.369466},   0.535778, {0.281271, 0.49604, 0.980923, 0.267273, 0.708723,    0.949374, 0.0306195, 0.746098, 0.392868, 0.677306, 0.425689,    0.917454, 0.203868, 0.134142, 0.373668, 0.498916, 0.903961,    0.390487, 0.0376524, 0.633368, 0.914573, 0.384288, 0.668186,    0.0975895, 0.633302, 0.888248, 0.687263, 0.830317, 0.924579,    0.938874}, {0.656644, 0.0842191, 0.531711, 0.261567, 0.230955,    0.166765, 0.327844, 0.127425, 0.857286, 0.667849, 0.423883,    0.736938, 0.819634, 0.034481, 0.50931, 0.35265, 0.151448,    0.936891, 0.876007, 0.464402, 0.464185, 0.106575, 0.951428,    0.525529, 0.807541, 0.0223557, 0.419717, 0.263961, 0.576586,    0.855591}, {0.0918732, 0.136536, 0.7193, 0.187742,    0.66799}, {0.399598, 0.899666, 0.153261, 0.158681, 0.0469482},   1, {{3, 3}, {2, 4}, {3, 4}, {3, 4}, {3, 0}, {1, 0}, {1, 1}, {3,     3}}, 0.21637,   0.282673}, {-22.672879543323493817, -0.62860637974940195936,{0.0000187062, 0.999629, 2.65313e-7, 0.000118007,    0.000233541}}}, {{{0.582546, 0.284033, 0.109795, 0.331245,    0.0570174, 0.476493, 0.0874391, 0.911528, 0.793056, 0.899907,    0.231848, 0.819655, 0.65652, 0.180607, 0.0441063, 0.151665,    0.256921, 0.280941, 0.890845, 0.992984, 0.209973, 0.532723,    0.674476, 0.710311},   0.627427, {0.24869, 0.564681, 0.379066, 0.57041, 0.772197,    0.477242, 0.467538, 0.777354, 0.872291, 0.245394, 0.647883,    0.120834, 0.691684, 0.201287, 0.496218, 0.863913, 0.410743,    0.310442, 0.503234, 0.65394, 0.87802, 0.635966, 0.792923,    0.0265126, 0.62933, 0.0712855, 0.413856, 0.456103, 0.857133,    0.594044}, {0.946318, 0.678749, 0.984842, 0.34865, 0.298435,    0.557915, 0.293158, 0.147363, 0.802216, 0.694002, 0.882415,    0.836921, 0.298982, 0.0400618, 0.00439587, 0.200955, 0.50606,    0.0135492, 0.375066, 0.129669, 0.0922033, 0.557447, 0.517934,    0.535625, 0.145885, 0.878698, 0.533092, 0.186975, 0.84745,    0.320783}, {0.239934, 0.0396124, 0.0452338, 0.626782,    0.357518}, {0.202691, 0.746251, 0.58672, 0.353122,    0.00173693}, -1, {{3, 2}, {2, 0}, {2, 0}, {3, 4}, {0, 1}, {1,     1}, {2, 2}, {3, 3}}, 0.573171,   0.978056}, {-14.903268081345429808, -1.2420724207377910018,{0.146834, 0.396591, 0.285168, 0.104488, 0.0669188}}}, {{{0.872068,    0.147988, 0.015724, 0.460123, 0.336442, 0.0021033, 0.137026,    0.927031, 0.149467, 0.154653, 0.816243, 0.687097, 0.109855,    0.109419, 0.189462, 0.329579, 0.907163, 0.363168, 0.602742,    0.976456, 0.905426, 0.122976, 0.0295712, 0.9984},   0.0333585, {0.974987, 0.0138472, 0.538277, 0.696916, 0.972884,    0.876821, 0.611246, 0.547449, 0.818231, 0.0605777, 0.924149,    0.437594, 0.708812, 0.871116, 0.59457, 0.530431, 0.345644,    0.268374, 0.618114, 0.625004, 0.222668, 0.238803, 0.619714,    0.591646, 0.247681, 0.224956, 0.0814369, 0.89473, 0.274797,    0.348135}, {0.470191, 0.347281, 0.456566, 0.287557, 0.546041,    0.909687, 0.747754, 0.416441, 0.951471, 0.379256, 0.402109,    0.148067, 0.333357, 0.754252, 0.179441, 0.909264, 0.713642,    0.162606, 0.93176, 0.684308, 0.632205, 0.267876, 0.656963,    0.336172, 0.162015, 0.920595, 0.200398, 0.048615, 0.615974,    0.0109085}, {0.452644, 0.632174, 0.664503, 0.631652,    0.0505347}, {0.484107, 0.331146, 0.877401, 0.871094, 0.574843},   0, {{2, 0}, {1, 3}, {1, 3}, {2, 0}, {2, 1}, {0, 0}, {3, 4}, {2,     2}}, 0.714795,   0.939334}, {-15.577470525768623627, -0.98541604207881939555,{0.999998, 3.22124e-7, 1.39541e-7, 9.33738e-7,    3.45103e-7}}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      NormalDistribution prior_logbeta = 
        normal_distribution (data[i].inputs.logbetamean, 
                             data[i].inputs.logbetastddev);

      float logbeta = data[i].inputs.logbeta;
      float pz[5];

      for (unsigned int j = 0; j < 5; ++j)
        {
          pz[j] = drand48 ();
        }

      float q = em (data[i].inputs.loga,
                    6,
                    &logbeta,
                    data[i].inputs.c,
                    5,
                    data[i].inputs.d,
                    pz,
                    data[i].inputs.logpriorz,
                    data[i].inputs.clamp,
                    data[i].inputs.ratings,
                    8,
                    &prior_logbeta.base,
                    1000,
                    1e-7);

      assert (fabs (q - data[i].desired_output.newq) <=
              1e-4 * (1 + fabs (q) + fabs (data[i].desired_output.newq)) ||
              ( fprintf (stderr, "%u: %g ?= %g\n", i, q, data[i].desired_output.newq),
                0 ));
      
      assert (fabs (logbeta - data[i].desired_output.logbeta) <=
              5e-4 * (1 + fabs (logbeta) + fabs (data[i].desired_output.logbeta)) ||
              ( fprintf (stderr, "(oops logbeta ... but this is ok) %g ?= %g\n", logbeta, data[i].desired_output.logbeta),
                1 ));

      for (unsigned int j = 0; j < 5; ++j)
        {
          assert (fabs (pz[j] - data[i].desired_output.pz[j]) <=
                  5e-4 * (1 + fabs (pz[j]) + fabs (data[i].desired_output.pz[j])) ||
                  ( fprintf (stderr, "(oops pz[%u] ... but this is ok) %g ?= %g\n", j, pz[j], data[i].desired_output.pz[j]),
                    1 ));
        }
    }
}

int 
main (void)
{
  srand48 (69);

  test_em ();

  return 0;
}
